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Abstract 

The consistency of second-order closure models with results from hydrodynamic stability 
theory is analyzed for the simplified case of homogeneous turbulence. In a recent study, 
Speziale, Gatski and MacGiolla Mhuiris [Phys. Fluids A 2 , 1678, 1990] showed that second- 
order closures are capable of yielding results that are consistent with hydrodynamic stability 
theory for the case of homogeneous shear flow in a rotating frame. It is demonstrated in 
this paper that this success is due to the fact that the stability boundaries for rotating 
homogeneous shear flow are not dependent on the details of the spatial structure of the 
disturbances. For those instances where they are - such as in the case of elliptical flows 
where the instability mechanism is more subtle - the results are not so favorable. The 
origins and extent of this modeling problem are examined in detail along with a possible 
resolution based on Rapid Distortion Theory (RDT) and its implications for turbulence 
modeling. 
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1. INTRODUCTION 


While Reynolds stress models are usually thought of as being completely empirical and 
without any theoretical foundation, recent work has demonstrated that this is far too pes- 
simistic an assessment of the current generation of models 1-3 . One encouraging result is the 
recent discovery that second-order closure models are capable of accurately predicting the 
stability boundaries for homogeneous shear flow in a rotating frame (see Speziale, Gatski and 
MacGiolla Mhuiris 4 ). This is a predictive capability that no eddy viscosity model possesses 
- including models ranging from the Baldwin-Lomax model 5 to the standard K — e model 6 . 
Such better predictions are usually thought to arise from the fact that second-order closures 
are based on the Reynolds stress transport equation which incorporates more turbulence 
physics since it is a rigorous consequence of the Navier-Stokes equations. However, it is 
important to understand to what extent these results genuinely have a sound theoretical 
basis. This forms the motivation for the present paper. 

It will be shown that the reason existing second-order closure models do so well for 
homogeneous shear flow in a rotating frame is due to the fact that the stability boundaries 
do not depend on the details of the spatial structure of the disturbances. For homogeneous 
turbulent flows where this is not the case - such as the example of elliptical flows which 
exhibit a more subtle instability mechanism involving resonance that only exists over a 
narrow band of wavenumbers - the results are not so favorable. This problem will be 
documented in the sections to follow along with a possible resolution. At least some limited 
consistency with results from hydrodynamic stability theory is needed if Reynolds stress 
models with greater predictive capabilities are to be obtained. Otherwise, it would not even 
be possible to predict whether a statistically unsteady turbulence decays or is self-sustaining. 
The level of consistency that Reynolds stress models should have with results from linear 
stability theory will be discussed in detail and some illustrative calculations will be presented. 

2. THEORETICAL BACKGROUND 

We will consider incompressible, homogeneous turbulent flows with constant mean veloc- 
ity gradients, dvi/dxj. The Reynolds stress tensor is a solution of the transport equation 7 


~/j — T~ik 
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for homogeneous turbulence in an arbitrary non-inertial reference frame where 0; is the 
angular velocity of the reference frame and is the permutation tensor. In (1), n^- and 
£ij are, respectively, the pressure-strain correlation and the dissipation rate tensor given by 
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where rq is the fluctuating velocity, p is the fluctuating pressure, v is the kinematic viscosity 
and the overbar represents an ensemble mean (t 8J = ufuf is the kinematic Reynolds stress 
tensor). Second-order closures in turbulence are based on the solution of a modeled version 
of (1) which arises from a second moment of the Navier- Stokes equations. 

In order to achieve closure, models for IRj and are needed. In virtually all existing 
second-order closures, these correlations are modeled in the form 


O — 

II ij = eAi 3 (h) + K Mi 3 ki(h)-^— (4) 
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where K = is the turbulent kinetic energy, e = is the (scalar) turbulent dissipation 
rate, and bij = (rp — ^K8ij)/2K is the Reynolds stress anisotropy tensor. Eq. (4) is based on 
a simplified analysis of the Poisson equation for the fluctuating pressure whereas (5) follows 
from the Kolmogorov assumption of local isotropy 7 . For two-dimensional mean turbulent 
flows that are homogeneous and near an asymptotic equilibrium state - which includes the 
cases to be considered in this paper - (4) simplifies to the form (see Speziale, Sarkar and 
Gatski 8 ) 
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are the mean rate of strain and absolute mean vorticity tensors and C\ — C 5 are constants. 

Virtually all existing second-order closures use a modeled transport equation for the 
turbulent dissipation rate that is of the general form 


. _ n £ dv t 
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for homogeneous turbulence. In (8), C e 1 and C e 2 are either constants or are functions of 
the flow invariants which become constants when the turbulence has achieved an asymptotic 
state. The Speziale, Sarkar and Gatski (SSG) model 8 , the Launder, Reece and Rodi model 9 
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and the IP model 10 are all special cases (6) and (8) where C e i = 1.44 and C e 2 is 1.83 for the 
former model and 1.92 for the latter two models. 


3. MODEL PREDICTIONS 


For homogeneous shear flow in a rotating frame (see Figure 1(a)) 


0 5 0 \ 

0 0 0 , Sl t = (0,0,0) (9) 

V 0 0 0 / 

where S is the shear rate and 0 is the angular velocity of the reference frame which are both 
constants. In homogeneous shear flows, the Reynolds stress anisotropy tensor bij and the 
dimensionless shear parameter SK/e achieve equilibrium values that are largely independent 
of the initial conditions (these fixed points are denoted by (&p)oo and (SK/e ) 0 0 since they 
are the asymptotic values obtained in the limit as t — > oo). From the contraction of (1), and 
the results of Speziale and MacGiolla Mhuiris 11 , it follows that for long times 



k* 


c e2 - c el 
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SK 
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after (8) is made use of (here, K* = K/ K 0 and K* is its time derivative with respect to 
the dimensionless time, t* = St). If (e/SK ) 0 0 > 0, then we have the long time asymptotic 
solution 


I<* ~ e At * (11) 

where 

A — C e 2 _ C el f £ ^ 

= c sl - 1 \~skJ 

is the dimensionless growth rate. This, of course, corresponds to unstable flow since there is 
an exponential growth of the turbulent kinetic energy given that C e i > 1 and C e2 > C e \. 

If (e/ SK ) oo = 0, then K* — > 0 as t — > oo in such a way that the long time asymptotic 
solution is the power law decay K* ~ t~ 13 where the exponent [3 depends on the ratio of 
the rotation rate to the shear rate, Si/ S, as well as on the model constants in (6) and (8) 
(see Speziale and MacGiolla Mhuiris 11,12 ). The fixed point (e/ SK ) 0 0 = 0 exists for all Si/ S, 
however it is only stable (i.e., a stable focus) outside of the interval 

A <^< B ( 13 ) 
where A and B are determined by the constants in the model. For A < SI/ S < B } a stable 
fixed point of the focus type, with (e/SK ) 0 0 > 0, exists in parallel with the (e/SK ) 0 0 = 0 
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fixed point that becomes a saddle which is unstable. This exchange of stabilities is illustrated 
in the bifurcation diagram provided in Figure 2. Virtually all existing second-order closures 
predict that point A is approximately at ft/S = —0.09. Furthermore, second-order closures 
with C 5 = 0 (such as the early proposal of Rotta 13 ) predict that point B is somewhat 
beyond ft/S = 0.5. The SSG model was optimized to perform well for a range of benchmark 
two-dimensional mean turbulent flows that are near equilibrium; this renders the constant 
C 5 = 0.40 which places point B at ft/ S ~ 0.53 (see Speziale, Sarkar and Gatski 8 ). However, 
it is worth noting that the older Launder, Reece and Rodi 9 model predicts a premature 
restabilization at ft/ S ~ 0.38. 

It is thus clear that state-of-the-art second-order closures predict that homogeneous shear 
flow in a rotating frame is unstable for approximately the range 

ft , , 

0 < - < 0.5 (14) 

O 

wherein there is an exponential growth of the turbulent kinetic energy. Discernibly outside 
of the interval (14) the flow is stable wherein the turbulent kinetic energy undergoes a 
power law decay (the time evolution of the turbulent kinetic energy in each of these distinct 
regimes is illustrated in Figure 3 for the SSG model). These results are consistent with linear 
hydrodynamic stability theory which predicts unstable flow for 0 < ft/S < 0.5 (see Lezius 
and Johnston 14 ). Of course, one expects the nonlinear stability boundaries to be somewhat 
broader as predicted by the models (it is well known that homogeneous shear flow, where 
ft/S = 0, is linearly stable but nonlinearly unstable 15 ). Unfortunately, a detailed comparison 
with results for the nonlinear stability of rotating homogeneous shear flow is not possible 
due to the fact that no such studies have apparently been published. Nonetheless, based on 
previous comparisons with large-eddy simulations 8 , the predictions of second-order closures 
appear to be very good, at least from a partial quantitative standpoint. 

While the stability boundaries of rotating homogeneous shear flow are described surpris- 
ingly well by second-order closures, the same is not the case for certain other homogeneous 
turbulent flows where the instability mechanism is more subtle. An example of such a ho- 
mogeneous turbulent flow is the case of plane strain with an added solid body rotation in an 
inertial frame (see Figure 1(b)). This is characterized by the mean velocity gradient tensor 


dvj 

dxj 


/ 0 T + lo 0 

T -uj 0 0 

V o oo 


(15) 


where T is the strain rate and u is the imposed rotation rate which are both constants (here, 
fti = 0 since the reference frame is inertial). For u > T, we have the case of elliptical flows 
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which have been the subject of recent studies in hydrodynamic stability theory (see Bayly 16 
and Landman and Saffman 17 ). We have hyperbolic flows for lo < T. 

For this flow, Eqs. ( 1 ) and (8) can be written in terms of K* and e/TK as follows: 


dK* 

dt* 


( 4^12 + 


(16) 
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where Ilf, = an d t* = Tt (here we do not need to write the transport equation for 

b 33 since b 33 = —611 — b 22 from the traceless property of the anisotropy tensor). Following 
the same analysis as presented in Speziale and MacGiolla Mhuiris 11 , the fixed points are 
obtained by setting the derivatives on the left-hand-sides of ( 17 )- ( 20 ) to zero and then solving 
the resulting algebraic equations simultaneously (here, the phase space is four dimensional). 
Neglecting the quadratic return term for simplicity, since it is small, yields: 

(^“ = -rG w)„ (21) 

and 
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. r k / 00 

for —6 < cu/r < 6 . Here, we have 
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given that a = (C e 2 — l)/((7 e i — 1); ( e/TK ) co = 0 outside of the interval —6 < u/T < 6. 
When (21) is substituted into (16) it yields an equation that is identical to (10) except that 
(e/SK ) 00 is replaced by ( e/TK ) QO which is given by (22)(the expression for (e/SK ) 0 0 valid 
in rotating homogeneous shear flow can be found in Speziale and MacGiolla Mhuiris 11 for 
comparison with (22)). 

Analogous to the case of homogeneous shear flow in a rotating frame, existing second- 
order closures predict unstable flow for the intermediate band of rotation rate to strain rate 
ratios, —6 < u /T < h, where ( e/TK ) QO > 0 yielding exponential growth (for the SSG model, 
6 ss 1.34). For u/T outside of this interval, ( e/TK ) QO = 0 yielding a power law decay of the 
turbulent kinetic energy. The fixed point ( e/TK ) QO = 0 exists for all values of u/T, but it 
is only stable outside of the interval —6 < u/T <6. This is illustrated by the bifurcation 
diagram (Figure 4) and the plots of the time evolution of the turbulent kinetic energy (Figure 
5) provided for the SSG model. There is the same type of exchange of stabilities (with an 
identical phase space portrait) as that obtained for the case of rotating homogeneous shear 
flow. While this stability picture is correct for homogeneous shear flow in a rotating frame, it 
is decidedly in error for plane strain with an added solid body rotation as given by (15). For 
\uj/T\ > 1 we have the class of elliptical flows which have been shown to be unstable for all 
finite values of u/T (see Bayly 16 , Landman and Saffman 17 , Pierrehumbert 18 and Waleffe 19 , 
the latter of which considers the nonlinear stability). All existing second-order closures 
erroneously predict a restabilization starting somewhere before |u;/r| = 2. The reason for 
these poor predictions in contrast to the remarkably good predictions for homogeneous shear 
flow in a rotating frame will be addressed in the next section. 

4. LINEAR STABILITY ANALYSIS 

A simplified linear stability analysis of the homogeneous flows considered in the previous 
section will now be presented. Our purpose here is to determine how sensitive the stability 
boundaries are to the details of the spatial structure of the disturbances. We consider a 
standard decomposition of the velocity held ry and pressure p of the form 

Vi = Ui + it-, p = P + p (24) 

where Ui and P represent the velocity and pressure associated with the base how whereas 
u\ and p' represent the disturbance velocity and pressure. The substitution of (24) into the 
Navier-Stokes equation, in a non-inertial frame, yields the disturbance equation 15 
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after the base flow (which also satisfies the Navier-Stokes equation) is subtracted out and 
terms that are nonlinear in the disturbances, which are assumed small, are neglected. 
Eq. (25) forms the basis for the linear theory of hydrodynamic stability in non-inertial 
frames. 

The simplified case of spatially uniform disturbances of the form 


>)=!■'(«), p=p(t) 


(26) 


will be considered. Here, u\ = u'(t) implies that p 1 = p'(t) after we make use of the in- 
compressibility constraint and the fact that the pressure disturbances must be bounded for 
the linear stability limit to be valid. After substituting (26) into (25), the simple dynamical 
system 

u t + A lk u k = 0 (27) 

is obtained where 


Atk -i : + 2e ' ik{1 ’ 


(28) 


is a constant tensor for the homogeneous flows under consideration. For homogeneous shear 
flow in a rotating frame, it is clear from (9) that 


/ 0 A -20 0 
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0 
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Since the non-zero eigenvalues of this matrix are given by 


A 1j2 = ±^20 (A - 20) 

it follows that the solution of (27) is of the form 


u ' a = A*(0)e 


Xrvt 


(30) 
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relative to the principal axes x a . Hence, the kinetic energy K of the disturbances has the 
asymptotic form 

(31) 

(32) 


where t* = St and 
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is the rotational Richardson number (see Bradshaw 20 who defined Ri to be the negative of 
(32)). Hence, with this definition of the Richardson number, the flow is unstable if Ri > 0 
or, equivalently, if 

0 1 

0 < 5 < 2 ; (33 > 

Eq. (33) is identical to the rigorous linear stability result based on a full normal mode 

analysis (see Lezius and Johnston 14 ). It is thus clear that here the instability mechanism 
does not depend on the details of the spatial structure of the disturbances - a result that 
explains the success of second-order closures which smear out information about the spatial 
structure of the fluctuations in homogeneous turbulent flows. 

For plane strain with an imposed solid body rotation in an inertial frame, the matrix 
takes the form 

/ o r + w o\ 

Ai k = T-lo 0 0 (34) 

V o oo/ 

since, for this case, = dUi/dxk ■ The non-zero eigenvalues of (34) are as follows: 

A 1j2 = ±Vr 2 -u 2 . (35) 

Thus, implementing the same logic used to arrive at (31), we obtain the asymptotic form 


K oc eVi-W 


where t* = Yt. Hence, we have unstable flow for |u;/r| < 1: a result that is comparable to 
the bound of |u;/r| < 1.34 obtained from the SSG second-order closure. However, unlike the 
corresponding result obtained for rotating homogeneous shear flow, this stability result is 
decidedly not correct. Recent linear stability studies based on a full normal mode analysis 
- conducted by Bayly 16 for the inviscid case and Landman and Saffman 17 for the viscous 
case - have indicated that there is unstable flow for all finite values of u/Y in the elliptical 
flow domain where |u;/r| > 1. Furthermore, it should be noted that this flow is nonlinearly 
unstable for all finite values of u/Y (see Waleffe 19 ). This is a somewhat surprising result 
in that strong rotation is often thought of as a stabilizing effect. However, for rotation 
dominated elliptical flows, the instability mechanism is based on a more subtle resonance 
effect that only exists over a narrow band of wavenumbers 16 . This feature causes considerable 
difficulties in conducting direct numerical simulations of elliptical flows (see Blaisdell and 
Shariff 21 ). 

Interestingly enough, second-order closures predict the same topological structure of the 
instability for plane strain in a rotating frame of reference as that manifested in the bifur- 
cation diagram shown in Figure 3 for the hyperbolic-elliptical flow case. However, unlike for 
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the latter case, this instability picture is correct for plane strain in a rotating frame which 
bears some similarities to rotating shear flow 22 . For plane strain in a rotating frame, 


m 

dxi 


/ o r o 

r o o | , n t = (o,o,o) 

V 0 0 0 


rendering 


/ 0 
= r + 20 
V o 

The non-zero eigenvalues of (37) are given by 


r - 20 o 
o o 
o o 


(37) 


a 1j2 = ±Vr 2 - 40 2 

which yields a real positive eigenvalue - and, hence, unstable flow - if |0/T| < |. This is in 
line with RDT results which predict a restabilization for values of |0/T| that exceed 0.7 or 
even slightly smaller values (F. Godeferd, private communication). The fixed point analysis 
for the SSG second-order closure in this case is qualitatively the same as that for plane strain 
with an added solid body rotation shown earlier. It yields unstable flow for 

0 

— 0.6 < — < 0.6 


which is in line with RDT and with the results of the simplified linear stability analysis 
presented above. The time evolution of the RDT solution for the turbulent kinetic energy is 
displayed in Figure 6 for O/T = 0.3 and O/T = 2 which, respectively, are in the unstable and 
stable regime (F. Godeferd, private communication). From a qualitative standpoint, these 
results show a remarkable similarity to the results obtained from second-order closures shown 
in Figure 5 for the companion problem of plane strain with an added solid body rotation. 
Second-order closures do not see the qualitative difference between plane strain in a rotating 
frame and plane strain with an added solid body rotation. 

Finally, we would like to discuss the role that Rapid Distortion Theory (RDT) can play 
in the development of Reynolds stress turbulence closures that are more consistent with 
hydrodynamic stability theory results for homogeneous flows. We first note that a general 
RDT analysis is equivalent to a linear stability analysis for homogeneous flows - a fact that 
is not often stated explicitly and, therefore, not widely appreciated outside of the held of 
theoretical turbulence. For homogeneous turbulence, 


On. 

dx, 


= 0 


(38) 
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and, thus, the evolution equation for the fluctuating velocity reduces to 

dui dui dui dvi dp 

~K7 + v kj, — = u k - U k - — 

OX OX}~ OX}- OX}- OXi 

in an arbitrary non-inertial frame. RDT deals with solutions to the linearized version of 
(39), where the nonlinear convective term u k dui/dx k is neglected 23 . As a result of these 
considerations, the Reynolds- averaged Navier-Stokes equation is identical to the Navier- 
Stokes equation itself eliminating any distinction between v k and U k . Under these conditions, 
(39) becomes identical to (25), which renders an RDT analysis equivalent to a linear stability 
analysis. 

It is thus clear that improved consistency with linear stability theory results would be 
achieved if Reynolds stress turbulence closures were more consistent with RDT. In homoge- 
neous turbulence, the RDT solution has been demonstrated to constitute an excellent ap- 
proximation to the full Navier-Stokes equations, for at least short elapsed times, in rapidly 
distorted flows where 

dvj .. £o_ 

dxj K 0 

given that || • || denotes any suitable norm. Unfortunately, within this domain where it is 
expected to be an excellent approximation, existing second-order closures have been shown 
to perform poorly compared to RDT. This has motivated Reynolds 24 and Kassinos and 
Reynolds 25 to attempt to develop Reynolds stress closures that are fully consistent with 
RDT by incorporating the structure tensor that includes information on the dimensionality 
of the turbulence. While this recent work of Reynolds and co-workers has a worthwhile 
goal - and contains many interesting new ideas - questions can be raised about whether 
or not it is possible to develop any one-point closure that is consistent with RDT for all 
flows. Two-point closures, which contain full spatial information, are probably needed for 
this purpose 26 . 

As stated earlier, it has been clearly demonstrated during the past decade that existing 
second-order closures perform poorly in the RDT limit. However, it must be noted that 
while the development of second-order closures that are consistent with RDT for the early 
transients of rapidly distorted turbulent flows would help, it is not sufficient to guarantee 
consistency with linear hydrodynamic stability theory. It is possible to construct a second- 
order closure that in the RDT limit has approximately the correct behavior but for weaker 
distortions behaves erroneously. Thus, the first author has recently pursued an approach 
that seeks to establish consistency with RDT only for a limited set of benchmark results, 
without compromising the performance of second-order closures in more weakly distorted 
turbulent flows where nonlinear effects play a major role 27 . 
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5. CONCLUSION 


An analysis of the consistency of existing second-order closure models with hydrodynamic 
stability theory results has been conducted for two basic homogeneous turbulent flows. Most 
notably in this regard, it is crucial for consistency to ensure that a turbulence model does not 
erroneously predict a decaying turbulent kinetic energy in a flow that is linearly unstable. In 
contrast to the earlier, and highly encouraging, results obtained for homogeneous shear flow 
in a rotating frame, it was found that the predictions are poor for the case of elliptical flows. 
The existing hierarchy of second-order closures predict that elliptical flows restabilize when 
the rotation rate exceeds the strain rate by any discernible amount whereas recent linear 
stability analyses indicate definitively that the flow is unstable for all finite rotation rates. 
It was found that this difficulty results from the fact that the instability mechanism is more 
subtle in elliptical flows. It encapsulates a resonance effect that exists only over a narrow 
band of wavenumbers which yields sensitivity to the spatial structure of the disturbances 

- information that is smeared out in homogeneous turbulence when Reynolds averaging is 
implemented. On the other hand, the stability boundaries for homogeneous shear flow and 
plane strain in a rotating frame are relatively insensitive to the spatial structure of the 
disturbances since they can be obtained from a stability analysis based on spatially uniform 
disturbances. 

It is argued that better consistency with the results of linear stability theory would be 
achieved if models were more consistent with RDT. However, full consistency with linear 
stability theory would require a model that is consistent with RDT but within a framework 
that yields growth rates for any unstable flow that are bounded below by the rapid distor- 
tion limit (milder distortions where nonlinear interactions play an important role should, in 
general, be less stable than the linear limit of RDT). As mentioned earlier, despite the very 
interesting work of Reynolds and co-workers 24,25 , it is debatable as to whether any one-point 
closure can be developed that is fully consistent with RDT for all flows. Nonetheless, it 
does appear to be useful to have at least some limited consistency with RDT in certain 
benchmark homogeneous turbulent flows (Speziale and Xu 27 have attempted a simpler ap- 
proach wherein a relaxation time and non-equilibrium Pade approximation is made around 
the equilibrium solution of the traditional hierarchy of second-order closures). This more 
limited approach may be more practical considering the fact that the inconsistencies are not 
that widespread. Homogeneous turbulent flows involving strains or shear in a rotating frame 

- as well as hyperbolic flows - are well described. A detailed comparison of a variety of 
second-order closures with DNS data bases for elliptical flows, hyperbolic flows and plane 
strain in a rotating frame will be the subject of a future paper. 
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LIST OF FIGURES 


Figure 1. Schematic of the flows considered, (a) Homogeneous shear flow in a rotating frame, 
and (b) Plane strain with an added solid body rotation. 


Figure 2. Bifurcation diagram for second-order closures corresponding to homogeneous shear 
flow in a rotating frame. 

Figure 3. Time evolution of the turbulent kinetic energy for homogeneous shear flow in a 
rotating frame: predictions of the SSG second-order closure model for an initial condition of 
SK 0 /e 0 = 4. 

Figure 4. Bifurcation diagram of the SSG second-order closure model corresponding to 
plane strain turbulence with an added solid body rotation (elliptical flows correspond to 

M/r| > i). 


Figure 5. Time evolution of the turbulent kinetic energy for plane strain turbulence with an 
added solid body rotation: predictions of the SSG second-order closure model for an initial 
condition of TK 0 /e 0 = 2. 

Figure 6. Time evolution of the turbulent kinetic energy for plane strain in a rotating frame: 
RDT calculations of F. Godeferd. 
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